home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / obsolete / anova_unequal.pro < prev    next >
Text File  |  1997-07-08  |  12KB  |  493 lines
  1. ; $Id: anova_unequal.pro,v 1.2 1997/01/15 04:02:19 ali Exp $
  2. ;
  3. ;  Copyright (c) 1991-1997, Research Systems Inc.  All rights
  4. ;  reserved. Unauthorized reproduction prohibited.
  5.  
  6.  
  7.  pro printout,TName,BName,IName,SST,SSB,SSI,SSE,R,C,N,FT,FB,FI,unit
  8.  
  9. if N_Elements( TName) EQ 0 THEN TName='Treatment'
  10. if N_Elements( BName) EQ 0 THEN BName='Block'
  11. if N_Elements( IName) EQ 0 THEN IName ='Interaction'
  12. printf,unit, format='
  13. printf,unit,'          Source        SUM OF SQUARES    DF     MEAN SQUARE       F       p'
  14.       printf,unit,'*******************************************************************************'
  15.  
  16. MSSE = SSE/(N-R*C)
  17. DFE = N-R*C
  18.  
  19. MT = SST/(C-1) 
  20. FT = MT/MSSE
  21. printf,unit,         $
  22.      Format='(A14,7X,G15.7,3X,I5,3X,G11.4,1X,G11.4,3X,G11.4)',$
  23. TName,SST,C-1,MT,FT, 1-F_Test1(FT,C-1,DFE)
  24.  
  25. MT = SSB/(R-1) 
  26. FB = MT/MSSE
  27. printf,unit,        $
  28.     Format='(A14,7X,G15.7,3X,I5,3X,G11.4,1X,G11.4,3X,G11.4)', $
  29.              BName,SSB,R-1,MT,FB, 1-F_Test1(FB,R-1,DFE)
  30.  
  31. MT = SSI/((C-1)*(R-1))
  32. FI = MT/MSSE
  33. printf,unit,         $
  34. Format='(A14,7X,G15.7,3X,I5,3X,G11.4,1X,G11.4,3X,G11.4)',  $
  35.        IName,SSI,(C-1)*(R-1),MT,FI, 1-F_Test1(FI,(C-1)*(R-1),DFE)
  36.  
  37. printf,unit,         $
  38.       Format='(A14,7X,G15.7,3X,I5,3X,G11.4)', 'Error',SSe,DFe,MSSe
  39.  
  40. RETURN
  41. END
  42.  
  43.  
  44. function mult_matrix ,X,Y,R,C
  45. ; Return the product of X with the design matrix
  46. ; specified by Y, Y(i,j) = the number of observations
  47. ; in cell (i,j). R and C are number of rows and columns
  48.  
  49. Here = 0
  50. First = X(0)
  51. X = X(1:*)
  52. Result = [0]
  53. skip = INDGEN(R-1)
  54. for i = 0L,R-1 DO  BEGIN
  55.  
  56.    for j = 0L, C-1 DO  BEGIN
  57.       
  58.         if j EQ C-1 THEN tot = -Total(X(R-1:C+R-3))      $
  59.          else tot = X(j+R-1)
  60.  
  61.        if i EQ R-1 THEN tot = tot - Total(X(0:R-2))          $
  62.           else tot = tot +  X(i) 
  63.        if j EQ C-1 and i ne R-1 THEN tot= tot -      $
  64.                       total(X(R+C-2 + i*(C-1): R+C-3+(i+1)*   $
  65.                            (C-1)))
  66.        if i EQ R-1 and j ne  c-1 THEN tot = tot -             $
  67.                                 Total(X(R+C-2 +j + skip*(C-1)))    
  68.        if(i EQ R-1 and j EQ C-1) THEN tot = tot +             $
  69.                                              Total(X(R+C-2:R*C-2))
  70.        if i NE r-1 and j NE C-1 THEN tot = tot + X(i*(C-1)+j+R+C-2)
  71.        
  72.        RESUlT = [RESULT, Replicate(tot,Y(i,j))]  
  73.    ENDFOR
  74.  ENDFOR
  75. X = [First,X]
  76. RETURN,Result(1:*) + First
  77. END
  78.  
  79.  
  80. Function Fit, G, R, C
  81.  
  82. ; Compute B such that X*B = G, where X is the design matrix
  83. ; for an experiment with R rows, C columns and one observation
  84. ; per cell
  85.  
  86. M = total(G)/(R*C)
  87. CF = Fltarr(R*C)
  88. CF(0) = M
  89. CF(1:R) = transpose((Replicate(1.0/C,C) # G) - M)
  90. CF(R:R+C-1) = G # Replicate(1.0/R, R) - M
  91. CF(R+C-1:*) = G(0:C-2,0:R-2) - replicate(1.0,C-1)#CF(1:R-1) - $
  92. CF(R:R+C-2)#replicate(1.0,R-1) - M
  93. return, CF
  94. END
  95.  
  96. function RSI_INV, G, R, C
  97. ; same as Fix but assume no interactions
  98. CF    = fltarr(R+C-1)
  99. CF(0) = Total(G)/(R*C)
  100. Row = Replicate(1.0,C) # G
  101. CF(1)     =  Total(Row(0) - Row)/(R*C)
  102.  
  103. if R gt 2 THEN CF(2:R-1 ) = CF(1) - (Row(0) - Row(1:R-2))/C
  104. Col =G #  Replicate(1.0,R)
  105. CF(R)     = Total(Col(0) - Col)/(R*C)
  106. if C gt 2 THEN $
  107.  CF(R+1:*) = CF(R) - (Col(0)-COl(1:C-2))/R 
  108.  
  109. return,CF
  110. END
  111.  
  112. function FitI,G,R,C
  113. ; regression fit to anova model with observations in G,
  114. ; no interactions and only 1 observation per cell.
  115. A = RSI_INV(G, R, C)
  116. CF = fltarr(R*C)
  117. CF(0:R+C-2) = A
  118. Y = fltarr(R,C) + 1
  119. CF = Mult_Matrix(CF,Y,R,C)
  120. return,CF
  121. END
  122.  
  123.  
  124. pro ComputeSSE, B ,YY,Y, R, C, I, SSE,BB
  125. ;ComputeSSE implements the generalized gradient algorithm
  126. ;in 'Nonorthogonal Analysis of Variance Using a Generalized
  127. ;Conjugate-Gradient Algorithm in Journ. Amer. Statistical
  128. ;Association, 1982, vol. 7.
  129.  
  130. BB   = FLTARR(R*C)             ;INITIALIZE
  131. G   = YY*Y
  132. SSE = Total(YY*G)
  133. RR    = FLTARR(R*C)
  134. eps  = 1.0
  135. d    = 0.0
  136. d1  = 1.0
  137. NUM = 0
  138. P = R * C
  139. while d1 GT 1e-8 DO BEGIN
  140.  NUM = NUM + 1
  141.  Z = FitI(G,R,C)
  142.  d1  = Total(Z*Z)      
  143.  eps = ABS(d1 - d)
  144.  if(NUM GT 1) THEN c1 = d1/d ELSE c1 = 0
  145.  d = d1
  146.  RR = Z + c1*RR
  147.  Z =  Y * RR
  148.  tt = Total(RR*Z)
  149.  if tt ne 0 THEN $
  150.    a = d/Total(RR*Z) $
  151.  else message, $
  152.   "Data causes division by zero. "
  153.  
  154.  BB = BB + a*RR
  155.  SSE = SSE - a*d
  156.  G = G - a * Z
  157. ENDWHILE
  158. BB = RSI_INV(reform(BB,c,r),r,c)
  159. return
  160. END
  161.  
  162. pro Make_Matrix,X,Y,R,C
  163. ; take anova data and make regression matrix X. Y(i,j) = number of observations
  164. ; in the ith,jth cell.
  165.  
  166.  
  167. X = fltarr(R*C,Total(Y))
  168. Here = 0
  169. X(0,*) = 1
  170. skip = INDGEN(R-1)
  171. for i = 0L,R-1 DO  BEGIN
  172.  
  173.    for j = 0L, C-1 DO  BEGIN
  174.  
  175.         TEMP = fltarr(R*C-1)    
  176.         if j EQ C-1 THEN temp(R-1:C+R-3) = -1 else temp(j+R-1) = 1
  177.        if i EQ R-1 THEN temp(0:R-2) = -1 else temp (i) = 1
  178.        if j EQ C-1 and i ne R-1 THEN          $
  179.                     temp(R+C-2 + i*(C-1): R+C-3+(i+1)*   $
  180.                                                     (C-1)) = -1 
  181.        if i EQ R-1 and j ne  c-1 THEN      $
  182.                           temp(R+C-2 +j + skip*(C-1))=-1    
  183.        if(i EQ R-1 and j EQ C-1) THEN temp(R+C-2:R*C-2) =1
  184.        if i NE r-1 and j NE C-1 THEN temp(i*(C-1)+j+R+C-2) =1
  185.        
  186.        X(1:*,HERE:HERE+Y(i,j)-1) = temp # Replicate(1.0,Y(i,j))  
  187.        HERE = HERE + Y(i,j) 
  188.    ENDFOR
  189.  ENDFOR
  190. RETURN
  191. END
  192.  
  193.  
  194.  
  195.  
  196.  
  197.  
  198.  
  199. pro make_matrix1,T,Y,R,C
  200. ;Compute X#transpose(X) directly
  201.  
  202. T = Fltarr(R*C,R*C)
  203.  
  204. LastRow = Total(Y(R-1,*))
  205. LastCol = Total(Y(*,C-1))
  206. T(0,0) =  Total(Y)            ;Multiply by column 0---all 1's
  207. for i = 1L, R*C-1 DO BEGIN
  208.  if (i LE R-1) THEN T(0,i) = Total(Y(i-1,*)) - LastROW    $
  209.  ELSE if (i LE R+C-2) THEN T(0,i) = Total( Y(*,i-R)) - LastCOL $
  210.  ELSE BEGIN
  211.   i1 = (i- R - C + 1)  
  212.   a  = i1 MOD (C-1)
  213.   b  = i1/(C-1)
  214.   T(0,i) = Y(b,a) - Y(R-1,a) - Y(b,C-1) + Y(R-1,C-1)
  215.  ENDELSE
  216. T(i,0) = T(0,i)
  217. ENDFOR
  218.  
  219.  
  220.  
  221. ;*************************************************************
  222. ;                  MULTIPLY BY BLOCKS
  223.  
  224. for i = 1L,R-1 DO  BEGIN
  225.  
  226.  for j= i,R-1 DO BEGIN        ;Block * BLOCK********
  227.   If (i EQ j) THEN T(i,i) = Total(Y(i-1,*)) + LastROW  $
  228.   ELSE T(i,j) = LastROW
  229.   T(j,i) = T(i,j)
  230.  ENDFOR
  231.  
  232.  for j = R,R+C-2 Do  BEGIN    ;Block * TREATMENT***********
  233.   T(i,j) = Y(i-1,j-R)- Y(R-1,j-R) -Y(i-1,C-1) + Y(R-1,C-1)
  234.   T(j,i) = T(i,j)
  235.  ENDFOR
  236.  
  237.  for j = R+C-1,R*C-1 Do BEGIN    ;Block * Interaction******
  238.   j1 = (j - R - C +1)
  239.   a  = j1 MOD (C-1)
  240.   b  = j1/(C-1)
  241.   if (b EQ i-1) THEN T(i,j) = Y(b,a) + Y(R-1,a) -   $
  242.                              Y(b,C-1) - Y(R-1,C-1)      $
  243.   ELSE T(i,j) = Y(R-1,a) - Y(R-1,C-1)
  244.   T(j,i) = T(i,j)
  245.  ENDFOR
  246.  
  247. ENDFOR
  248.  
  249. ;**************************************************************
  250. ;        MULTIPLY BY BLOCKS
  251.  
  252. for i = R,R+C-2 Do BEGIN           ;TREATMENT * TREATMENT ***************
  253.     i1 = i - R
  254.  for j = R,R+C-2 Do BEGIN
  255.   j1 = j - R
  256.   if (i EQ j) THEN T(i,j) = Total(Y(*,i1)) + LastCOL  $
  257.   ELSE T(i,j) = LastCOL
  258.   T(j,i) = T(i,j)    
  259.  ENDFOR
  260.  
  261.  for j = R+C-1,R*C-1 Do BEGIN      ;BlOCK*INTERACTION
  262.   j1 = (j - R - C +1)
  263.   a = j1 MOD (C-1)
  264.   b = j1/(C-1) 
  265.   if ( a EQ i1) THEN T(i,j) = Y(b,a) - Y(R-1,a) + Y(b,C-1)   $
  266.                              - Y(R-1,C-1)                   $
  267.   ELSE T(i,j) = Y(b,C-1) - Y(R-1, C-1)
  268.   T(j,i) = T(i,j)
  269.  ENDFOR
  270. ENDFOR
  271.  
  272. ;**************************************************************
  273. ;              MULTIPLY BY INTERACTION
  274.  
  275. for i = R + C -1, R*C -1 DO BEGIN
  276.  i1 = (i - R - C + 1)
  277.   a1 = i1 MOD (C-1)
  278.   b1 = i1/(C-1)
  279.  for j = i, R*C-1 DO BEGIN
  280.   j1 = (j - R - C + 1)
  281.   a = j1 MOD (C-1)
  282.   b = j1/(C-1)
  283.  if (i EQ j) THEN T(i,j) = Y(b,a) + Y(R-1,a) + Y(b,C-1) +  $
  284.                            Y(R-1,C-1)                      $
  285.  ELSE if ( a EQ a1) THEN T(i,j) = Y(R-1,a) + Y(R-1,C-1)    $
  286.  ELSE if ( b EQ b1) THEN T(i,j) = Y(b,C-1) + Y(R-1,C-1)    $
  287.  ELSE T(i,j) = Y(R-1,C-1)
  288.  T(j,i) = T(i,j)
  289.  ENDFOR
  290. ENDFOR
  291.   
  292.     
  293. RETURN
  294. END 
  295.  
  296.  
  297. function mult,B,Y,R,C,I
  298. ; Return the product of X with Transpose( design matrix
  299. ; specified by Y, Y(i,j) = the number of observations
  300. ; in cell (i,j)). R and C are number of rows and columns
  301. Tot = Total(Y)
  302. Tot1 = Total(Y(R-1,*))
  303. LastRow = Total(B(Tot-Tot1 : *))
  304. B1 = [Total(B)]
  305. Next = -1 
  306.  
  307. for i = 0L, R-2 DO BEGIN
  308.  Last = Next + 1
  309.  Next = Next + Total(Y(i,*))
  310.  Temp = Total( B(Last:Next)) - LastRow
  311.  B1 = [B1,Temp]
  312. ENDFOR
  313.  
  314. Tot = 0
  315. TotUp = Fltarr(R,C)
  316. for i = 0L,R-1 DO       $
  317.  for j = 0L,C-1 DO BEGIN
  318.   TotUp(i,j) = Tot
  319.   Tot = Tot + Y(i,j)
  320.  ENDFOR
  321.  
  322.  
  323. TEMP1 = [0]
  324. for j = 0L,C-1 DO BEGIN
  325.  Temp =0
  326.  for i = 0L, R-1 DO                    $
  327.   TEMP = TEMP + Total(B(TotUP(i,j):TotUp(i,j) + Y(i,j) - 1))
  328.  if ( j NE C-1) THEN                  $
  329.   TEMP1 = [TEMP1,TEMP]                
  330. ENDFOR
  331.  
  332. TEMP1 = TEMP1(1:*)
  333. TEMP1 = TEMP1 - TEmp
  334. B1 = [B1,TEMP1]
  335.  
  336. for i = 0L,R-2 DO                     $
  337.  for j = 0L, C-2 DO BEGIN
  338.   TEMP = Total(B(TotUP(i,j): TotUp(i,j) + Y(i,j) -1))
  339.   TEMP = TEMP - Total(B(TotUp(R-1,j):TotUp(R-1,j) +  Y(R-1,j)-1))
  340.   TEMP = TEMP - Total(B(TotUP(i,C-1):TotUp(i,C-1) + Y(i,C-1)-1))
  341.   TEMP = TEMP + Total(B(TotUp(R-1, C-1): TotUp(R-1,C-1) +      $
  342.               Y(R-1,C-1)-1))
  343.   B1 = [B1, TEMP]
  344.  ENDFOR
  345.           
  346.  
  347.  
  348. Return,B1
  349. END
  350.           
  351.  
  352.  
  353.  
  354.  
  355. pro anova_unequal, Data1,FT_Test,FB_Test,FI_Test, TName = TN, $
  356.                   BName = BN, IName = IN,Missing = M,    $
  357.                   ListName = LN, No_Printout =NP
  358. ;+
  359. ;
  360. ; NAME:
  361. ;    ANOVA_UNEQUAL
  362. ;
  363. ; PURPOSE:
  364. ;    Perform two way analysis of variance with interactions and unequal
  365. ;    cell sizes - that is, not every treatment / block combination
  366. ;    has the same number of iterations.
  367. ;
  368. ; CATEGORY:
  369. ;    Statistics.
  370. ;
  371. ; CALLING SEQUENCE:
  372. ;    ANOVA_UNEQUAL, Data, FT_Test, FB_Test, FI_Test
  373. ;
  374. ; INPUTS: 
  375. ;    Data:    Array of experimental data, dimensioned (Treatment#, I, B),
  376. ;        where I is the number of repetitions in each cell.
  377. ;
  378. ; OUTPUT:
  379. ;    Anova table displaying  Sum of Squares, Mean Squares, and F Ratio 
  380. ;    for sources of variation.
  381. ;
  382. ; OPTIONAL OUTPUT PARAMETERS: 
  383. ;      FC_Test:    value of F for treatment or column variation.
  384. ;
  385. ;      FB_Test:    value of F for row or block variation.
  386. ;
  387. ;      FI_Test:    value of F for column variation.
  388. ;
  389. ; KEYWORDS:
  390. ;      MISSING:    missing data value.  If undefined, assume no missing data. If 
  391. ;        unequal sample sizes, set M to place holding value.
  392. ;
  393. ;    LIST_NAME:    name of output file. Default is to the screen.
  394. ;
  395. ;    TNAME:    name to be used in the output for treatment type.
  396. ;
  397. ;    BNAME:    name to be used in the output for block type.
  398. ;
  399. ;    INAME:    name to be used in the output for interaction type
  400. ; 
  401. ;  NO_PRINTOUT:    flag , when set,  to suppress printing
  402. ;                       of output
  403. ; RESTRICTIONS:
  404. ;    Each treatment/block combination must have at least 1 observation.
  405. ;
  406. ; SIDE EFFECTS:
  407. ;    None.
  408. ;
  409. ; PROCEDURE:
  410. ;    Overparameterized least squares procedure with sum to zero 
  411. ;    restrictions.
  412. ;- 
  413.  
  414. On_Error,2
  415.  
  416.  
  417. if N_Elements(LN) NE 0 THEN openw,unit,/get,LN else unit = -1
  418. if N_elements(data1) eq 0 THEN BEGIN
  419.    printf,unit,"anova_unequal- Data array is empty"
  420.    return
  421. ENDIF
  422.  
  423. SD = size(Data1)
  424. R  = SD(3)
  425. C  = SD(1)
  426. IN  = SD(2)
  427.  
  428. if SD(0) ne 3 THEN BEGIN
  429.    printf,unit, "anova_unequal- Data array must be 3 dimensional."
  430.    return
  431. ENDIF
  432.  
  433. DATA =DATA1
  434. B = [0]
  435. Y = Fltarr(R,C) + R
  436. YY = Y
  437.  
  438. for i = 0L,R-1 DO         $
  439.  for j = 0L,C-1 do BEGIN
  440.   if( N_Elements(M)) THEN BEGIN
  441.    here = where(Data(j,*,i) NE M, count)
  442.    Y(i,j) = count
  443.  
  444.    if count EQ 0 THEN BEGIN
  445.     printf,unit,'anova_unbalanced - stop,empty cell'
  446.     goto, done
  447.    ENDIF
  448.  
  449.    here = transpose(Data(j,here,i))
  450.    B = [B,here]
  451.    YY(i,j) = Total(here)
  452.    ENDIF ELSE BEGIN
  453.     here = transpose(Data(j,*,i))
  454.     B = [B,here]
  455.     Y(i,j) = IN
  456.     YY(i,j) = Total(here)
  457.    ENDELSE
  458.  ENDFOR 
  459.  
  460. if N_Elements(M) EQ 0 THEN  $
  461.    printf,unit,'anova_unequal, there is no missing data'
  462.  
  463. YY = YY * 1.0/Y
  464. B = B(1:*)
  465.  
  466. BB = fit(transpose(YY),R,C)
  467. SSE = Total((B-mult_matrix(BB,Y,R,C))^2)
  468.  
  469. W =   1.0/Y # replicate(1.0,C)
  470. W = 1.0/(W/(C^2))
  471. YYY =   YY # replicate(1.0/C,C)
  472. WA = Total(w*YYY)/total(w)
  473. SSB = Total(w*(YYY-WA)^2)
  474.  
  475. W =    replicate(1.0,R) # (1.0/Y)
  476. W = 1.0/(W/(R^2))
  477. YYY = replicate(1.0/R,R) # YY
  478. WA = Total(w*YYY)/total(w)
  479. SST = Total(w*(YYY-WA)^2)
  480.  
  481. ComputeSSE, B ,transpose(YY),transpose(Y), R, C, I, SSI,BB
  482.  
  483. if( N_Elements(NP) EQ 0) THEN     $
  484. printout, TName, BName, IName, SST,SSB,SSI,SSE,R,C, Total(Y),  $
  485.  FT_Test,FB_Test,FI_Test,unit
  486.  
  487. DONE:
  488. if (unit NE -1) THEN Free_Lun,unit
  489.  
  490. return
  491. END
  492.  
  493.